    %%simulation sec4
function y = getTraj()
    global N M In Ifo IMN r cord%swarm para
    global L L1 Ld L1d Lup %topolohy
    global LR L1R LupR %repaired topolohy
    global A B A0 B0 Bbar%state matrix para
    global T K1 %virtual layer para
    global invPR PR Zetadot Zdot X0dot fol_a%private mat Epshat Epshatdot
    global Kcp Hcp P %CP layer param
    global igx export tgx%middle var
    global DoSP DoSR rho rho_R phi rhogx%DoS attacks
    global dbar kesen % threshold of Actuation atatcks
    global gbar1 tau_Ra
    global LupA LupRA ttg1x ttg1y tridx1
    global Zetaall tau_Zbar
    global eZ1 eZ2 eZ3 EpsZ1 EpsZ2 EpsZ3 XiZkeep h
    global EpsZ1bar EpsZ2bar EpsZ3bar gg1x gg2x gg3x tritime1x tritime2x tritime3x fuzhichi
    global tkp1ag1x tkag1x tkp1ag1y tkag1y tkp1ag2x tkag2x tkp1ag2y tkag2y tkp1ag3x tkag3x tkp1ag3y tkag3y
    global fismat %Fuzzy policy

    %read config files
    % 获取当前脚本所在目录，即项目目录下MatlabSrc/C4-main0526
    coder.extrinsic('fileparts')
    coder.extrinsic('fullfile')
    coder.extrinsic('readtable')
    currentDir = fileparts(mfilename('fullpath'));
    % 构建配置文件所在目录路径
    configDir = fullfile(currentDir, '../../Config');

    % 读取attack_strategy_config.csv文件
    attackFilePath = fullfile(configDir, 'attack_strategy_config.csv');
    attackTable = readtable(attackFilePath);

    % 提取并处理attack_strategy_config.csv中的数据
    % dos_period
    dosPeriodStr = attackTable.dos_period(1);
    dosPeriodTokens = strsplit(char(dosPeriodStr), ',');A
    dosPeriod = zeros(6, 2);
    for i = 1:length(dosPeriodTokens)
        rangeTokens = strsplit(dosPeriodTokens{i}, '-');
        dosPeriod(i, 1) = str2double(rangeTokens{1});
        dosPeriod(i, 2) = str2double(rangeTokens{2});
    end
    % dos_freq
    dosFreq = attackTable.dos_freq(1);
    % impersonation_signal
    impersonationSignal = attackTable.impersonation_signal(1);
    % fdi_signal
    fdiSignal = attackTable.fdi_signal(1);
    % aa_signal
    aaSignal = attackTable.aa_signal(1);
    % aa_rate
    aaRate = attackTable.aa_rate(1);

    % 读取defense_strategy_config.csv文件
    defenseFilePath = fullfile(configDir, 'defense_strategy_config.csv');
    defenseTable = readtable(defenseFilePath);

    % 提取并处理defense_strategy_config.csv中的数据
    % event_trigger_slope
    eventTriggerSlope = defenseTable.event_trigger_slope(1);
    % zeno_free_interval
    zenoFreeInterval = defenseTable.zeno_free_interval(1);
    % backup_topology_time
    backupTopologyTime = defenseTable.backup_topology_time(1);
    % controller_gain_k
    controllerGainKStr = defenseTable.controller_gain_k(1);
    controllerGainKStr = controllerGainKStr{1};
    controllerGainKTokens = strsplit(controllerGainKStr, ',');
    controllerGainK = zeros(1, 2);
    controllerGainK(1) = str2double(controllerGainKTokens(1));
    controllerGainK(2) = str2double(controllerGainKTokens(2));
    % controller_gain_h
    controllerGainHStr = defenseTable.controller_gain_h(1);
    controllerGainHStr = controllerGainHStr{1};
    controllerGainHTokens = strsplit(controllerGainHStr, ',');
    controllerGainH = zeros(1, 2);
    controllerGainH(1) = str2double(controllerGainHTokens(1));
    controllerGainH(2) = str2double(controllerGainHTokens(2));


    %fismat=readfis('4');
    % fuzhichi=0.01*sqrt(5);
    fuzhichi=aaSignal*sqrt(5);
    tkp1ag1x=0; %tkag1x=0;
    tkp1ag1y=0; tkag1y=0;
    tkp1ag2x=0; tkag2x=0;
    tkp1ag2y=0; tkag2y=0;
    tkp1ag3x=0; tkag3x=0;
    tkp1ag3y=0; tkag3y=0;
    tritime1x=[];
    tritime2x=[];
    tritime3x=[];
    gg1x=0;
    gg2x=0;
    gg3x=0;
    tridx1=0;
    ttg1x=[];
    ttg1y=[];
    eZ1=[]; eZ2=[]; eZ3=[]; EpsZ1=[]; EpsZ2=[]; EpsZ3=[];XiZkeep=[];
    EpsZ1bar=[]; EpsZ2bar=[]; EpsZ3bar=[];

    rho=0; rho_R=0;
    kesen = 0.1;
    dbar = 2;
    igx = 1;
    Zetadot = zeros(12,1);
    Zdot=Zetadot;
    X0dot = zeros(4,1);
    % Epshatdot = zeros(12,1);
    %hall = zeros(16,1);
    rhogx = [];
    export = [];
    tgx = [];
    r = 2;
    cord = 2;
    DoSP = dosPeriod;

    % 5.5 7;
    %     8 9.5;
    In=eye(cord);

    M = 1;
    N = 3;
    Ifo=eye(N);
    IMN=eye(M+N);

    L=[0 0 0 0;
    -1 1 0 0;
    0 -1 1 0;
    0 0 -1 1;];

    scrB=-L(M+1:end,1);
    L1=L(M+1:end,M+1:end);
    Lup = L(M+1:end,:);%
    LupA=[0 0 0 0;0 0 0 0;0 0 0 0];
    %
    ups=inv(L1')*ones(N,1);
    Theta=diag(ups);
    invTheta = inv(Theta);
    Omega = (Theta*L1+L1'*Theta);

    kappa0 =eigs(inv(2*Omega)*diag(ups),1);
    %tau2 = 0.5;

    %
    L1d=kron(L1, eye(cord));
    Ld=kron(L, eye(cord));
    %

    LR = [0 0 0 0;
    -1 1 0 0;
    0 -1 1 0;
    0 -1 0 1];
    L1R=LR(M+1:end,M+1:end);
    LupR = LR(M+1:end,:);%
    LupRA =[0 0 0 0;0 0 0 0;0 0 0 0];
    upsR=inv(L1R')*ones(N,1);
    ThetaR=diag(upsR);
    invThetaR = inv(ThetaR);
    OmegaR = (ThetaR*L1R+L1R'*ThetaR);

    kappaR =eigs(inv(2*OmegaR)*diag(upsR),1);
    kappa = max(kappa0,kappaR);%*(gamma12)
    %pause
    %%

    A0=[0 1;-.5 0];%eig:+i,-i
    eig(A0)
    b=1;
    B0=eye(2)*b;
    A=[0 1;-1 -1.2];
    B=[0;1];
    Bbar = [1 0];
    PRIa=1;%杩涓や釜板互浠婚?,绉
    PRIb=0.4;%杩涓や釜板互浠婚?
    PRIc=PRIb*A0(2,1);
    PRId=PRIa+PRIb*A0(2,2);


    PR =[PRIa PRIb;
        PRIc PRId];%rank(PRI)=2蹇椤绘弧绉
    %PR = eye(2);
    invPR = inv(PR);

    %virtual layer parameter
    rho=0;% attack flag锛?0跺绘ㄧ
    phi = 0; % repair flag, 0跺绘ㄤ慨?T=eye(2);

    %tau_a = 1/0.4;
    % Ff=0.1;
    Ff=dosFreq;

    gamma1 =  fdiSignal;%浜0.15/sgx
    sgx =  0.5;


    syms gvar
    Asym=sym(A0);
    Bsym=sym(B0);
    Psym=plya(Asym,Bsym,gvar);%sqrt(gamma12)*

    Pkap = double(subs(Psym,gvar,gamma1));

    gamma1tilde = gamma1*sgx;%>0.09Ff=0.1
    gamma1bar = gamma1 - gamma1tilde;
    gamma2 =0.94; %LMIsolve

    flaggamma2 = (eig(A0*Pkap+Pkap*A0'-gamma2*Pkap)<0);%1涓gamma2

    kappa1 = 2*norm(A0,2);
    kappa2 = norm(kappa*(B0'*B0*Pkap),2);


    gamma1star = gamma1bar*min(eigs(kron(ThetaR,(Pkap)),1,'SM'),eigs(kron(Theta,(Pkap)),1,'SM'));%
    %gamma1star = min(eigs(kron(Theta,(Pkap)),1,'SM'),eigs(kron(ThetaR,(Pkap)),1,'SM'))

    theta1mini = gamma1star*min(eigs((inv(L1))'*inv(L1),1,'SM'),eigs((inv(L1R))'*inv(L1R),1,'SM'));%

    Gm = Pkap;%inv(Pkap)
    theta2max = 2*kappa*max(eigs(kron(Omega,b^2*Gm*Gm),1),eigs(kron(OmegaR,Gm*Gm),1))-gamma1star;
    % theta3max = max(norm(kron(inv(L1)'*Theta*(L1),kappa*Gm*Gm)-gamma1star*kron(inv(L1)',eye(2)),2),...
    %                norm(kron(inv(L1R)'*ThetaR*(L1R),kappa*Gm*Gm)-gamma1star*kron(inv(L1R)',eye(2)),2));
    GX1 = kron(inv(L1)'*Theta*(L1),kappa*b^2*Gm*Gm)-gamma1star*kron(inv(L1)',eye(2));
    GX2 = kron(inv(L1R)'*ThetaR*(L1R),kappa*b^2*Gm*Gm)-gamma1star*kron(inv(L1R)',eye(2));
    theta3max = max(eigs((GX1+GX1')/2,1),...
                eigs((GX2+GX2')/2,1));
    %c1 = 1.1*theta3max/theta1mini
    %c1star=1600;
    c1star = (theta3max^2+sqrt(theta3max^4+theta1mini*theta2max*theta3max^2))/(theta1mini*theta3max);
    %c1star =292;%%%solve by copying theta1mini theta2max theta3max to U1113.m
    flagc1star = (c1star>theta3max/theta1mini);%1涓gamma2
    gbar =sqrt((theta1mini-theta3max/c1star)/(theta2max+c1star*theta3max));
    % gbar1 = gbar*0.6
    % gbar2 = gbar*0.8;
    gbar1 = gbar*eventTriggerSlope;
    gbar2 = gbar*(sqrt(1-eventTriggerSlope^2));
    %gbar2 = 0.05
    % tau_Zbar = 1/kappa1*log(kappa1/kappa2*sqrt(gbar2)/(sqrt(M))+1)
    tau_Zbar = zenoFreeInterval;
    mu_1 = eigs(ThetaR*inv(Theta),1);
    mu_2 = eigs(Theta*inv(ThetaR),1);
    %mu = max(mu_1,mu_2);
    mu = (mu_1*mu_2)/2;
    tau_Zuse=floor(tau_Zbar*100)/100;
    % tau_R = (gamma1tilde - Ff*log(2*mu))/((gamma1tilde+gamma2)*Ff)-(tau_Zuse)
    tau_R = backupTopologyTime;
    K1 = kappa*b*Pkap;
    tau_Ra = floor(tau_R*10)/10; %<tau_R?
    DoSR = DoSP+kron(tau_Ra*ones(size(DoSP,1),1),[1 0]);%

    %%
    save('DOS','DoSP','DoSR')

    %CP layer parameter
    Qre=eye(2);
    U=1;
    [P,~,~]=care(A,B,Qre,U);
    %Bbar=[1 0];
    Gammabar=A0-A;
    Gamma = Gammabar(2,:);

    % Kcp=-U*B'*P;
    % Hcp=Gamma-Kcp;

    Kcp = controllerGainK;
    Hcp = controllerGainH;

    X0(1:2)=[-1 0.5];
    X0(3:4)=[0.3 0.6];
    X0(5:6)=[0.9, 0.9];
    X0(7:8)=[0 0];
    X0(9:10)=[-0.2, -0.8]*1;
    X0(11:12)=[0 0];
    X0(13:14)=[0.3, -0.7]*2;
    X0(15:16)=[0 0];%16*1
    xhat0=[1,1,0.3,0.8,  0.4,0.6,0.6,0.4,  0.8,0.5,0.5,0.5];%12*1

    rhohat0 = 0.1*ones(N*cord, 1);
    Z10=[X0(1:4) xhat0]';%16*1

    Z20=[X0(1:4) xhat0]';%16*1
    Z30=[X0(1:4) xhat0]';%16*1
    X1=[xhat0 X0 rhohat0' Z10' Z20' Z30']'*10;%涓茶;X0'

    Zetaall = [X0(1:4) xhat0]';

    tic;
    h=0.0004;
    a=0;
    b = 50;
    tspan = a:h:b;
    [t, y] = runge_kutta2(@myssFun0526, X1, h, a, b); % 忽略 t，仅保留 y

%    % 获取当前脚本所在目录
%    currentDir = fileparts(mfilename('fullpath'));
%    % 定义 CSV 文件名
%    csvFilename = fullfile(currentDir, 'trajectory.csv');
%
%    % 直接保存 y 矩阵（要求 y 是二维数值型矩阵）
%    csvwrite(csvFilename, y);
end